function ivp=rk4s(a,c,Za,h,fname)

% M?todo de Runge-Kutta para sistemas de ecuaciones 
% 
% Entrada
%         - a , b extremos del intervalo 
%         - Za=[x1(a), ..., xn(a)] condiciones iniciales
%         - M es el n?mero de subintervalo en (a,b)
% Salida  - T el vector de valores equiespaciados
%		de la variable independiente
%         - Z=[x1(t) . . . xn(t)] donde xk(t) es la 
%	    aproximacion a la soluci?n para la k_esima
%  		funci?n inc?gnita.
% La funci?n f(t,x) deber? ser ingresada en la l?nea de comando
% de Octave.

T=a:h:c
M=length(T);
Z=zeros(M,length(Za));
Z(1,:)=Za;
global b;
global I;
global K;
%KArr = [0:0.1:1];
K = 1;
%K = KArr(1);
b = 0.02;
I = 0.04;

for j=1:M-1   
	k1=h*feval(fname,T(j),Z(j,:));   
	k2=h*feval(fname,T(j)+h/2,Z(j,:)+k1/2);   
	k3=h*feval(fname,T(j)+h/2,Z(j,:)+k2/2);   
	k4=h*feval(fname,T(j)+h,Z(j,:)+k3);   

	Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
end

ivp = [T' Z]